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Abstract. The density matrix renormalization group (dmrg) is applied to some one-dimensional reaction- 
diffusion models in the vicinity of and at their critical point. The stochastic time evolution for these 
models is given in terms of a non-symmetric "quantum Hamiltonian" , which is diagonalized using the 
dmrg method for open chains of moderate lengths (up to about 60 sites). The numerical diagonalization 
methods for non-symmetric matrices are reviewed. Different choices for an appropriate density matrix in the 
non-symmetric dmrg are discussed. Accurate estimates of the steady-state critical points and exponents 
can then be found from finite-size scaling through standard finite-lattice extrapolation methods. This is 
exemplified by studying the leading relaxation time and the density profiles of diffusion-annihilation and 
of a branching- fusing model in the directed percolation universality class. 
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1 Introduction 

The density matrix renormalisation group (dmrg) tech- 
nique was invented by White jj] in 1992 as a new tool for 
the diagonalization of quantum chain spin Hamiltonians. 
The dmrg allows to study much larger systems than is 
possible with standard exact diagonalization methods and 
provides data with remarkable accuracy. The method has 
since been applied to a large variety of systems in quan- 
tum || and classical g physics. These applications usu- 
ally considered hermitian quantum Hamiltonians for spin 
chains or symmetric transfer matrices for two-dimensional 
classical systems. In both cases the dmrg is known to work 
very well. For a collection of reviews see jl. 

Motivated by the success obtained for these problems, 
some efforts have been recently devoted to problems in- 
volving diagonalization of non-symmetric matrices which 
are often encountered in various fields of physics; as ex- 
ample we mention the cases of low temperature thermo- 
dynamics of spin chains (q| or out-of-equilibrium classical 
systems ||[7|. 

The models investigated in this paper belong to the lat- 
ter class of problems. We consider one-dimensional reaction- 
diffusion systems where the stochastic time evolution of 
the system is given by a non-hermitean operator H to 
which we apply the DMRG algorithm. One of the mod- 
els considered displays a non-equilibrium phase transition 
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which is expected to be in the universality class of directed 
percolation, a typical class for many out-of-equilibrium 
systems and for which the critical exponents are known 
numerically to a high degree of accuracy. 

Our main motivation for this study is as follows. At 
equilibrium, studies of finite systems through diagonal- 
ization of a transfer matrix and/or a quantum Hamilto- 
nian and analysis of the results through finite-size scal- 
ing, eventually combined with precise extrapolation algo- 
rithms, has been one of the standard methods of studying 
phase transitions in systems with many strongly interact- 
ing degrees of freedom |pPjlC|]. Compared to more stan- 
dard diagonalization methods, the DMRG does allow to 
study much larger systems, even if it works best in situa- 
tions when the ground state of H is separated by a finite 
gap from the first excited state, that is, at a finite dis- 
tance away from the critical point. That advantage is at 
least partially offset by the fact that for reasons of nu- 
merical accuracy, it is preferable to apply the dmrg to 
models with open boundary conditions, which lead to es- 
timates of the critical parameters which converge usually 
slower than those found for periodic boundary conditions. 
Here we study to what extent a finite-size diagonaliza- 
tion study of non- equilibrium systems is feasible and in 
particular, to what degree of precision critical points and 
critical exponents can be estimated. For that purpose, we 
shall consider some critical non-equilibrium models with 
known properties to compare with our dmrg data. 
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Diagonalizing non-symmetric matrices is numerically 
much more demanding than for symmetric matrices. Po- 
tentially, numerical instabilities might arise at several stages 
of the calculation. We shall discuss in some detail various 
diagonalization methods for non-symmetric matrices and 
have tested our results throughout by working with two 
different diagonalization algorithms. 

We shall not only consider the calculation of eigen- 
values, which determine the relaxation times, but also of 
eigenvectors, which are needed in the calculation of matrix 
elements, as they arise for example in the calculation of 
density profiles. In general, we find that for non-symmetric 
matrices the numerical accuracy is not as good as for sym- 
metric ones. However, we shall show that the dmrg is well 
capable to accurately determine the values of the critical 
parameters. This makes this technique a useful general- 
purpose method for the study of non-equilibrium critical 
phenomena. 

This paper is organized as follows: in section |2j, we in- 
troduce reaction-diffusion systems and their description 
in terms of non-symmetric evolution matrices, in section 
H we review the basics of the DMRG algorithm, and present 
in section || a brief review of diagonalization methods 
for non-symmetric matrices. More specific details about 
choosing a convenient density matrix for non-equilibrium 
systems are described in section a. In sections O and u\ 
we present our numerical results concerning the critical 
properties of the model and then conclude. 



2 Reaction-diffusion processes 

We consider a chain of length L in which each site is either 
occupied by a particle (A) or empty (0). The time evolu- 
tion of the system is given in terms of microscopic rules 
involving only a pair of particles on neighbouring sites. We 
shall consider the following system, using the notation of 



A0 <-> %A (with rate D) 

AA -» 00 (with rate 2a) 

AA -> $A, A% (with rate 7) 

AH), <DA -> 00 (with rate S) 

A%, $A -> A A (with rate /3) 



(1) 
(2) 
(3) 
(4) 
(5) 



, all the other pro- 

) or an increase (ra) 

for a list of common 



Apart from the diffusion process 
cesses involve either a decrease 
in the number of particles (see [D I 
alternative notations). 

Once the reaction rates are given, the stochastic evolu- 
tion follows from the master equation, which can be writ- 
ten as 

d\P(t)) 



dt 



-H\P{t)) 



(6) 



where \P(t)) is the state vector. The elements of the quan- 
tum "Hamiltonian" H are given by 



(a\H\r) 

(a\H\a) 



-w(t — > a) ; a ^ t 



^w(cr^r) 



(7) 



where \a), \t) are the state vectors of the particle config- 
urations a, t and w(t — > a) denotes the transition prob- 
ability between the two states and is easily constructed 
from the rates (0^) of the elementary processes. Since H 
is a stochastic matrix, the left ground state (s\ is 



E< 



(8) 



with ground state energy Eq = 0, since (s\H — 0. All 
other eigenvalues £$ of H have a non-negative real part 
$lEi > Jl2j . Since the formal solution of eq. (ph is 



\P(t)) 



,-Ht 



P(0)> 



(9) 



the system evolves towards its steady state |P(oo)). Let 
r := infj 5R£"i for i 7^ 0. Often, one simply has 



P — E\ — Eo — Ei 



(10) 



since Eq = 0. If r > 0, the approach towards the steady 
state is characterized by a finite relaxation time r = 1/T, 
but if r = 0, that approach is algebraic. This situation is 
quite analogous to non-critical phases (r 7^ 0) and critical 
points (t = 0), respectively, which may arise in equilib- 
rium quantum Hamiltonians. 

It is clear that H is non-symmetric if w(a — > r) 7^ 
w(t — > a). However, if the detailed-balance condition 



w{a -► t)P s ({t}) = w{t -> <r)P s (M) 



(11) 



is satisfied, where P s ({a}) is defined by 
|P(oo)) = J^ (7 P s ({cr})|tT), His similar to a non-stochastic, 
but symmetric matrix K, without affecting the locality 
of the interactions, e.g. ]1^ , ^3{ . Detailed balance always 
holds when besides diffusion only a single reversible reac- 
tion is present. For several reversible reactions, the cases 
when detailed balance holds are given in |L4|, eq. (4.48). 
Only if detailed balance holds, the right ground state \s) = 
|P(oo)) is related to the known left ground state (s\ in a 
simple way |]l5| . 

In this paper, we shall study systems without detailed 
balance, since it is our aim to explore the DMRG in a 
setting as different from equilibrium physics as possible. 
Specifically, we shall consider two models: 

(1) Diffusion-annihilation. The rates are 



D = 2a=p , (5 = 7 = 5 = 



(12) 



This model is characterized by an algebraic approach to- 
wards the steady state with a gap vanishing in the ther- 
modynamic limit as T ~ l/L 2 . It offers to us the addi- 
tional advantage that H can be diagonalised through free- 
fermion methods (e.g. p4| ) and thus analytical results are 
available for comparison with the numerical data. 
(2) Branching- fusing. Here we take 



D = 2a = "f = 5 = I — p , f3 = p 



(13) 



This model presents a non-equilibrium phase transition. 
For p small, the annihilating processes dominate and the 
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Table 1. Some critical exponents for ID directed percolation, 
as obtained from series data Jig] and Monte Carlo simulations 
| [17|jl8| . In the last column, the values obtained by density ma- 
trix renormalization in the present work are listed. The num- 
bers in brackets give the estimated uncertainty. 



that H as defined on a open chain with L sites and has 
the local structure 



exponent 


series 


simulation 


this work 





0.27647(10) 


0.27649(4) 




0/v± 


0.2520(1) 


0.25208(4) 


0.249(3) 


Pi 


0.7338(1) 


0.728(1) 




0t/v± 


0.6690(1) 


0.664(7) 


0.667(2) 


"II 


1.7338(1) 


1.73383(3) 




Ki 


1.0969(1) 


1.09684(1) 


1.08(2) 


Q = v\\/v±. 


1.5806(2) 


1.58074(4) 


1.580(1) 



steady state is simply the empty lattice. However, if p 
becomes sufficiently large, the steady state contains par- 
ticles at some finite density. The transition between these 
two phases is expected to fall into the directed percola- 
tion universality class (see also appendix [A|). Although 
even in ID, there is no analytical information available, 
series expansion studies and Monte Carlo simulation have 
over the years yielded extreme ly p recise estimates of crit- 
ical exponents, as reviewed in M] an d |jj],[l8|]- We collect 
some of their results in table [j]7but skip over repeating 
the precise exponent definitions here. We use the notation 
9 = v\\/v±_ for the dynamical exponent in order to avoid 
confusion with the exponent z = 2/9 which is also often 
used. If D = 0, the model would reduce to standard site- 
bond directed percolation, but the presence of D should 
not change the universality class. The last column of table 
jy shows the values of the exponents obtained from finite- 
size scaling extrapolation of our numerical data. 

In both models (||Jl|), the empty lattice |0 ... 0) is ob- 
viously an adsorbing state and consequently also a steady 
state. Since H is real, its eigenvalues are either real or 
occur in complex conjugate pairs. The first excited eigen- 
value Ei is always real. 

For the branching-fusing model, we also considered the 
case of particle injection at the boundaries of the system, 
i.e. we added the reaction 







— ► A (with rate p') 



(14) 



at the two sites at the edges. In this situation, |0 . . . 0) is 
no longer a stationary state of the system since the lattice 
is occupied by a non-vanishing density of particles for all 
values of p. We analyze the shape of the density profiles 
as function of the parameter p and of the injection rate p' . 



3 Density matrix renormalization group 
algorithm 

We recall briefly in this section the basic structure of the 
DMRG algorithm JD. The task is to find selected approxi- 
mate eigenvalues and eigenvectors of a given Hamiltonian 
H. That desired eigenvector \tj;) is called a target state and 
the process of selecting \ip) is referred to as targeting (it is 
possible to target several states, see below). One assumes 



H 



L-l 



1,4+1 



(15) 



where /ij,i+i is a local Hamiltonian acting on a pair of 
nearest-neighbour sites (this condition of locality can be 
somewhat relaxed but we restrict here to the simplest 
case) . In this paper, we always consider free (open) bound- 
ary conditions, as typically done in dmrg studies. 

The DMRG is an iterative method: it produces a cho- 
sen eigenvector (usually, one targets the ground state or 
the first excited state) and its eigenvalue, starting from 
the target state of a small chain which is known from 
exact diagonalization of the Hamiltonian (|l5|), and then 
using it to find \ip) on chains with an increasing number of 
sites. This is made possible by projecting at each iteration 
step the full vector space to a smaller space where only 
a selected number m of states is kept. This projection is 
carried out via the density matrix as described below. 

A DMRG calculation, at least in the context of appli- 
cations to critical phenomena, proceeds in two steps. The 
first one is the infinite system method (ism) which we now 
describe. Suppose we are interested in the ground state of 
H. As the starting point, consider a chain of four lattice 
sites which can be represented as B 



(i). 



Bi. , where • de- 



>(i) 



notes a single site and -ET / are blocks at the left and right 
side of the chain. Initially, they contain only one spin, that 
is B r } = • (of course, the calculation may be started at 
larger lattices) 

At this point, the main loop begins. The Hamiltonian 
H is easily written down and its ground state wave func- 
tion ipo(&hihJriPr) can be found via standard diagonal- 
ization routines, where ai and j3 r denote degrees of free- 
dom of the blocks B r ' and B\ and the indices ii,j r refer 
to the spin degrees of freedom of the single lattice points 
in the middle of the chain. The density matrix for the left 
part of the system is defined as 

p {l) {a U ii;^ h ki) = ^ i>o(®l,ihjr, /3r)lpo(lhh,jr, Pr) 

(16) 
which we shall write in a short-hand notation as p = 
tr (|^o)(V'o|)) where tr denotes a partial trace either in 
the left or right part of the system (we shall reconsider 
the precise choice of p in section g). Next, one solves the 
eigenvalue problem p\fii) = u)i\Qi). The eigenvalues of 
the density matrix are non-negative and can be ordered 
according to u)\ > 0J2 > ^3 > . . . > 0. Furthermore, if 
the ground state vector of H is normalized according to 
(ipo\ipo) = 1) one has Y^i^i = 1. Each eigenvalue u>i is 
equal to the probability of finding the left part of the chain 
in the corresponding density matrix eigenvector |J?j) when 
the whole system is in the ground state V'o)- The config- 
urational space reduction is obtained by keeping only the 
first m dominant density matrix eigenvectors \fii) with 
i = 1, 2, ...m, corresponding to the m largest Wj. For- 
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mally, the truncation can be represented by 



O 



«'•) 



O rt 



B 



(2) 



In these steps, on the B r are updated according to eq. (170, 
while the Bi are taken from the blocks calculated previ- 
ze ously. Finally, the Bi are updated again through the se- 



where O m — [|J?i), . . . , |J? TO )]. The accuracy of the projec- 
tion operation can be described by the truncation error: 



e = 1 - ^2 w« 



(18) 



The projection operation is repeated for the right part as 
well to obtain B r ' (if there is a left-right symmetry, as 

(2) 

in the models studied here, B r ' is obtained simply from 

(2) 

reflection of B\ ). Performing these calculations, sparse- 
matrix diagonalization techniques will be needed to obtain 
the ground eigenstate |Vo) 0I H. On the other hand, since 
the whole spectrum of the non-sparse matrix p is required, 
this is best obtained via some standard routines. 

Combining the two blocks with new sites one gets 



,(2) 



,(2) 



B\ • »B r , e.g. a chain of L = 6 sites after the first 
pass through the main loop. The next pass through the 
main loop begins by writing down H for this longer chain. 
Applying this procedure repeatedly at the left and 
right part of the system, one generates larger systems. At 
each iteration step, two new sites are added in the middle 
of the chain and the boundaries are pushed further away 
from each other. Schematically, this may be illustrated as 



R (1) ..B, 



(i) 



BP..B 



2) 



B 



(L/2-1) 
I 



• •B 



(L/2-1) 



The ISM procedure is repeated, typically until L w 1000, 
and if there is a finite gap in the low-lying spectrum of 
H, finite-size effects can then be neglected. However, it is 
well known that the ISM method alone is not enough to 
yield precise results in the thermodynamic limit L — > oo 
for systems close to criticality |19],[I| . 

In such cases (e.g. pQ]), a second step in the dmrg 
calculation is required. It is best to use the finite system 
method (fsm) designed by White [jjj to accurately deter- 
mine properties of systems of finite length. The starting 
point of the fsm is the target vector \ip) for a chain of 
given length L, as generated by the ISM described above. 
At this point further iterations are started. First, one cal- 
culates better approximations for the blocks on the left 
part rep resenting more than L/2 — 1 sites, using as before 
eq. (|17]), while for the blocks on the right part, one uses 
blocks generated in previous iterations in order to keep 
the total length of the system fixed at L. Schematically 
this looks as follows 



B[ L/2 - l) . .B^ 2 ^ 



B[ L/2) . .Bi L ' 2 -^ 



B 



(L-3) 
I 



• •B™ 



Second, this procedure is reversed and the larger blocks 
on right part of the system are refined. Schematically, 



B\ L - 3) . .BP - B[ L - 2) . ,B^ 



B, (1) . .Bi L -^ 



quence 



B\ 1] . .B^ - B [2) . .Bi L ^ 



. . . - B\ L/2 - l) . .Bi L ' 2 -V 

until one is back at the left-right symmetric partition. The 
target vector extracted at this stage can be used as start- 
ing point for the next FSM iteration. 

These "sweeps" improve the results both for eigenval- 
ues and eigenvectors m. For a given lattice size, two or 
three sweeps are usually enough to achieve convergence. 

In practice, to calculate critical exponents, one needs 
data for chains of various lengths in order to perform a 
finite-size scaling analysis. For better efficiency, we used 
repeatedly the fsm to calculate from the same run quan- 
tities for chains of different lengths as follows. The blocks 
generated at the end of the fsm for a chain of length Lq 
are used as starting point for further DMRG calculation: 
first, we use the ism to enlarge the system is symmetri- 
cally until a length L\ > Lq is reached. Second, the FSM 
sweeps are started again. In this way we produced accu- 
rate results for systems of various lengths L < L\ < . . ., 
typically equally spaced, in the same dmrg run. 



4 Diagonalization methods for non-symmetric 
matrices 

In "conventional" DMRG calculations the computationally 
dominant part is the determination of (ground state) wave 
functions given by the eigenvectors of a large sparse sym- 
metric N x N matrix A using a diagonalization algorithm 
for large sparse symmetric matrices such as the Lanczos 
algorithm. Using the fundamental operation 



\w)=A\v), 



(19) 



one builds iteratively small tridiagonal n x n matrices 
T n = QnAQ n with n <?C N, where Q n is a N x n matrix 
with orthonormal columns. The matrices T n are diagonal- 
isable by standard techniques; it can be shown that their 
extreme eigenvalues converge for rather small n to the ex- 
treme eigenvalues of A. Moreover, their eigenvalues form 
Sturm chains: if A^ are the eigenvalues of T n in ascending 
order and (x% those of T n+1 , one has fXi < A, < Hi+i- This 
ensures monotonous, easily controlled convergence of the 
extreme eigenvalues. Q n relates the eigenvectors of T n to 
those of A: |Ao)a = Qk|Ao)t„- In the case of large sparse 
non-symmetric matrices, as they arise in transfer matrix 
dmrg and non-hermitian dmrg, the situation is much less 
satisfying, as the lack of symmetry leads to intrinsic prob- 
lems of numerical stability. 

Essentially, there are two algorithms available, the Ar- 
noldi algorithm and the non-symmetric Lanczos algorithm 
[ pi| . In our calculations, we used both methods, to check 
them against each other to guard against numerical fal- 
lacies. Usually, they were in excellent agreement, though 
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the professional package available for the Arnoldi algo- 
rithm p2| usually seemed to be somewhat more accurate 
than our self-made implementation of the non-symmetric 
Lanczos algorithm. 

Both also build on the fundamental operation 



the tridiagonal matrix: 



and its transpose 



\w) = A\v) 



(v\A 



(20) 



(21) 



in the case of the non-symmetric Lanczos algorithm. Again, 
one builds iteratively small n x n matrices with n <C N di- 
agonisable by standard techniques for small non-symmetric 
matrices. In both cases, one finds that the extreme eigen- 
values of the small matrices do not converge to those of 
the big matrices in the very systematic manner of the 
symmetric case, i.e. they do not form Sturm chains. 

The Arnoldi method iteratively generates one sequence 
of orthonormal vectors \qi) forming the columns of a ma- 
trix Q n = [|gi), . • • , \q n )] and a Hessenberg matrix H n = 
Q^AQ n with elements hij. One starts from a random vec- 
tor |<7o) of unit length and h\o — 1 and iterates using 



lk+i,k 



\q k+1 ) = A\q k ) 



k 



km 



(22) 



with hij = (qi\A\qj). h k+ i tk is determined by enforcing 
unit length for \q k+ i). Here, (%| = \qi) T ■ 

The non-symmetric Hessenberg matrix is then diag- 
onalized using the standard QR algorithm. Some of the 
eigenvalues of H n will converge against some of those of 
A, and associated eigenvectors of H n can be transformed 
into those of A using Q n . However, it is quite intricate to 
assure that eigenvalue convergence actually happens nu- 
merically. We will not discuss this issue further, as there 
is a freely available highly sophisticated Arnoldi package 
(arpack) [|H . 

Let us now discuss the non-symmetric Lanczos algo- 
rithm, which is quite easily implemented, but has to be 
carefully refined for numerical stability. Then it provides 
a highly satisfying alternative approach portable to ma- 
chines where arpack might not be available. 

Let us first summarize the algorithm as presented in 
Ref. plj . The non-symmetric Lanczos algorithm generates 
iteratively two sequences of vectors from whom eventually 
left and right eigenvectors are built. In an adaptation of 
the symmetric case, one forms matrices from these vec- 
tors, Q n = [\qi),... ,\q n )] and P n = [(pi|,. . ., (p n \], to 
generate an (incomplete) basis transformation from the 
N x N matrix A to a tridiagonal n x n matrix T„, 



•Ln ^nf^cr. 



a\ 71 . 
01 a 2 72 







0n 







7n-l 
-1 Oi n 



(23) 



and demands biorthogonality, P n Q n — I n - One then ob- 
tains the following equations determining the elements of 



7fc(Pfe- 



\t k ) = {A- a k I)\q k ) - 7fe-i|gfc-i) (24) 
(s k \ = (Pk\(A - a k I) - /? fc _i<p fc _i| (25) 



with 



a k = (p k \A\q k ) 
0k = \\\t k )\\ 

Ik = (sk\t k )//3k- 



(26) 

(27) 
(28) 



Note that the choice of the off-diagonal coefficients is not 
unique. One starts with \qo) = and (po| = and two 
non-orthogonal random vectors \qi) and {pi\, normed such 
that |gi) has unit length and (j)i\qi) — 1 as input for the 
first iteration. This leads to an iteratively growing T n . The 
iterations are continued, until the lowest (two) eigenvalues 
of T n are sufficiently converged. From the left and right 
eigenvectors of T n one forms the left and right eigenvectors 
of A: \X) A = Q„|A) T „ and {\\ A = {X\ Tn P n - 

Note that the non-symmetric Lanczos algorithm yields 
left and right eigenvectors on an equal footing, while the 
Arnoldi algorithm has to be run twice for A and A T (how- 
ever, with less time-consuming matrix-vector multiplica- 
tions). 

In the present form, the non-symmetric Lanczos al- 
gorithm is capable of yielding rather good eigenvalues. 
However, eigenvectors are determined with rather poor 
accuracy only. This leads to a poor choice for the density 
matrix, in turn to a non-optimized decimation and then, 
after some DMRG steps, to a noticeable degradation of the 
eigenvalues too. There are two reasons for this. 

(1) The construction of the two Lanczos vector se- 
quences guarantees them to be biorthogonal, or (pi\qj) = 
6ij . While true mathematically, this biorthogonality is nu- 
merically only true locally, i.e. if i and j are close. Globally, 
small overlaps develop, in particular when the extreme 
eigenvalues of T n start to converge. This loss of orthog- 
onality introduces numerical errors in the mapping from 
the eigenvectors of T n to those of A. An analogous loss of 
orthogonality occurs in the Arnoldi algorithm. 

(2) Even the symmetric Lanczos algorithm shows the 
phenomenon that if one wants a non-extreme, say the sec- 
ond, eigenvalue to converge, it may occur that the eigen- 
value converges to its true value, but if one enforces too 
strict convergence criteria, "jumps" and will converge to 
the lowest eigenvalue, too, see pG| for an example. This 
phenomenon is much more frequent in the non-symmetric 
Lanczos algorithm and because of the non-monotonic con- 
vergence behaviour much harder to detect and handle nu- 
merically. One therefore has to relax convergence criteria, 
which lead to a not so good eigenvalue, but in particular 
a rather bad eigenvector. 

Both issues can be addressed in a satisfactory way. Ba- 
sically, all Lanczos methods are clever implementations of 
the power method, which obtains the largest eigenvector of 
a matrix by applying A successively to a random start vec- 
tor \v): for n -> oc, \X max ) = A"|u)/||yl n |«)||. The disad- 
vantage of the power method is its very slow convergence; 
however, it accumulates no numerical inaccuracies from 
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previous iterations: it essentially restarts at each iteration 
with a better guess for the eigenvector. A power iteration 
is thus a suitable eigenvector "beautifier": once Lanczos 
has effectively converged, it can be used to improve the 
generated eigenvector (and also somewhat the eigenvalue) 
by eliminating accumulated numerical inaccuracies. As we 
try to improve the eigenvectors of the smallest eigenvalues, 
we apply the power method to A — > — A+ (tr A)I, turning 
the smallest eigenvalues into the largest by absolute value, 
as the real parts of all eigenvalues are non-negative, while 
conserving the eigenvectors. As we have a left and right 
eigenvector on equal footing, one has to slightly modify 
the power method as follows: Let |ro) and (l$\ be the right 
and left eigenvectors obtained by Lanczos for the small- 
est (now largest) eigenvalue. One now forms successively 
\r) = A | ro), A = {lo\r) / {Iq\tq) as improved eigenvalue, and 
sets \r±) = \r) / X as improved right eigenvector. Now one 
forms (l\ = (lo\A, X = (l\ri) / (l \ri) as improved eigen- 
value, and sets (li\ = (l\/X as improved left eigenvector, 
and restarts with the new eigenvector pair, until the eigen- 
value is satisfactorily converged. We find empirically that 
if we apply this procedure to the second eigenvalue in the 
spectrum there is no problem that it will start to converge 
to the first one for a small number of iterations. It is there- 
fore suitable to improve both the lowest (largest) and the 
second eigenvalue and the associated eigenvectors. 



5 The choice of the density matrix 

An important question to address in the case of a DMRG 
calculation for a non-symmetric problem concerns the choice 
of density matrix. For reference, we recall the situation 
found in studies of equilibrium quantum spin chains at 
a non-zero temperature. There, a non-symmetric transfer 
matrix is generated by applying a Trotter decomposition 
along the imaginary time direction. The density matrix 
typically used in these systems is defined by || 



,.(0 



/.M\ 



where \ip\ ) and \ip\ ) denote the (normalized) left and 
right eigenvectors corresponding to the i-th eigenvalue of 
the non-symmetric Hamiltonian defined in (jm . This choice 
is easy to implement numerically, since one avoids all prob- 
lems related to the possibility of complex eigenvalues of a 
non-symmetric density matrix, but has also a more pro- 
found justification. Recall that in White's argument jjj, 
the choice of density matrix does not rely on the Hamilto- 
nian being symmetric. For a symmetric Hamiltonian the 
density matrix obtained from the combination |-0o)(' ! / ; o| a l~ 
lows the construction of a trial ground state function |^o) 
whose distance from the exact ground state |^o) is min- 
imal p[. In the non-symmetric case, the density matrix 
defined by eq. (pQ) provides a basis set which minimizes 
simultaneously the distance of the trial vectors from the 

exact right and left eigenstates \ip\ ) and \tpf ), see ap- 
pendix pi 

We also tried out some alternatives for a symmetric 
density matrix. For example, one might consider p being 
defined by using the right eigenstate only 



tr 



{\4 ) ){4 ) \} 



(31) 



(this density matrix was used in a study of the g-symmetric 
Heisenberg chain [[7| ) . Alternatively, one can use a density 
nratrix with mixed terms 



J3] .. 



a,- tr 



(l^ + l^^l + ^l)], (32) 



with on a normalization constant. The merit of this density 
nratrix would be that, while keeping the advantages of 
being symmetric, it contains terms \ip\ /{ipi I which holds 
the relevant information in the case of DMRG quantum 
thermodynamics. 

Which of the density matrices pM @, pW © and 
p[ 3 l (32) is the 'best' one, depends on the Hamiltonian 



tr 



K><4 r) i} 



(29) 



where tr denotes the partial trace on part of the system, 
the superscripts I and r label the left and right eigenvec- 
tors of the transfer operator T (throughout this section, 

the line vector (?/>o | is the transpose of the column vector 

\iJjq ) and so on). The choice ( p9J ) is justified as the one 
that maximizes the partition function which describes the 
thermodynamics of the system. Since p is non-symmetric, 
it may have complex eigenvalues. In practice, however, 
non-real eigenvalues appear only after a certain number 
of dmrg steps. They are thought to come from numerical 
round-off errors. Several approaches have been invented 
to avoid these problems [0 . 

In the present study, we restrict ourselves to symmetric 
density matrices. Usually, convergence was best for the 
type 



under study and also of the properties of the eigenstates 
one wants to calculate. This choice is a priori a rather 
difficult one. Our aim is to select the density matrix that 
produces the most stable numerical results. How to do this 
is best illustrated by some examples. 

Table || shows the finite system method (fsm) itera- 
tions for the branching-fusing model defined in eq. (|13| ) 
for a chain of length L = 20, with m — 32 states kept and 
with p — 0.84, thus in the vicinity of the critical point (we 
estimate p c = 0.84036(1) in the next section). The data 
were obtained by using only the first excited state as a 
target state. The L = 20 lattice was constructed through 
ISM iterations which are not shown. Although L = 20 is 
a rather small value for dmrg calculations the data in 
table |2| illustrate a typical behaviour which is found for 
larger systems as well. 

The first two columns in table |] show the lengths Li 
and L r of the right and left parts, respectively, used to 
form the system of total length L = 20. The third and 
forth columns show the value of the gap r := E\{p,L), 

obtained from the density matrices p\ and p\ , as de- 



,M 



tr 



{l^WfVl^Wl} (30) 



fined in (|30|) and (|31|). In both cases the iterations dur- 
ing the fsm improve the estimation of the gap as it can 
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Table 2. Gaps r = E\(jp,L) for the branching-fusing model 
eq. ([13|) with p — 0.84 and m = 32 states kept as found from 
the iterations of the finite system method using the density 



matrices p\ ' 



and pj 



[2] 



The total length of the chain is kept 
fixed to Li + L r = 20, where Li and L r indicate the lengths of 
the left and right parts, which vary during the application of 
the finite system method. 



L, 


L r 


r from p[ 


r from p\ 


10 


10 


0.0109823881100 


0.0109738583392 


11 


9 


0.0109794248185 


0.0109738587466 


12 


8 


0.0109784192015 


0.0109738588520 


13 


7 


0.0109781170856 


0.0109738588725 


14 


6 


0.0109781024661 


0.0109738588779 


13 


7 


0.0109781050213 


0.0109738588777 


12 


8 


0.0109781103648 


0.0109738588780 


11 


9 


0.0109781100799 


0.0109738588781 


10 


10 


0.0109738162569 


0.0109738594171 


11 


11 


0.0109738097605 


0.0109738594170 


12 


12 


0.0109738065720 


0.0109738594168 


13 


7 


0.0109738036253 


0.0109738594168 


14 


6 


0.0109738015566 


0.0109738594168 


13 


7 


0.0109738011287 


0.0109738594168 


12 


8 


0.0109738036758 


0.0109738594169 


11 


9 


0.0109738061907 


0.0109738594168 


10 


10 


0.0109737989313 


0.0109738594171 



model defined in (EL3|) for chains of lengths up to L = 18, 
with m = 32 states kept and with p = 0.8403578, very 
close to the the critical point. Here, both ground and 
first excited state were targeted at the same time, us- 
ing both the "mixed" density matrix p^ and the un- 
mixed density matrix pW (with the shorthand notation 

p^ :— (p Q + p[ )/2). To estimate the effect of numer- 
ical inaccuracies in the diagonalization routines, the cal- 
culations using pW were redone using extended precision 
arithmetic, with 30 instead of 14 mantissa digits. That is 
clearly enough to preclude the possibility of arithmetic er- 
rors to the quoted numerical precision. We observe that 
the choice of the unmixed density matrix is clearly supe- 
rior, unlike the case of transfer matrix DMRG. 

Table [| shows the finite system methodjFSM) itera- 
tions for the diffusion-annihilation model (fL2[) for chains 
with up to L = 18 sites, with m — 32 states kept and p — 
1. Data were obtained by targeting the third eigenstate 
(lowest excited state), using both p^ and pM, and are 
compared with the exact result (see the following section) . 
The choice of the symmetric density matrix is clearly supe- 
rior over against the unsymmetric choice, by several digits 
of precision, even for the small number of dmrg iterations 
performed; the difference in precision increases even fur- 
ther for longer systems. 



be seen from the general trend of convergence. Also, it 
is apparent that consideration of the symmetric partition 
Li = L r = L/2 = 10 achieves a particularly large increase 
in precision. In practice, this partition should be used for 
the final estimates from the dmrg algorithm. Both esti- 
mates of the gap are quite close to each other, however 
the density matrix ( |30| ) provides the more stable results, 
with a convergence up to the 11th digit. 

Although the density matrices are constructed only 
from the first excited state as a target state, we found that 
also the estimates of the ground state energy are rather 
accurate i.e. Eo(p,L) « 10~ 8 for the density matrix p\ ' 

and Eo(p, L) w 10~ 4 for p\ . Both values are close to the 
exact result E (p, L) — 0, and show again that the density 

matrix p\ not only shows better convergence for the tar- 
geted first excited state, but also a much more accurate 
value for Eo(p, L) which was not targeted. This is a practi- 
cally important remark for the following reason. Suppose 
one wants to find the first excited state. In cases with a 
very small gap r, there might be a spurious interchange 
with the ground state if the the error of the DMRG is larger 
than the true value of r. 

Targeting simultanenously the ground and first ex- 
cited states, i.e. using density matrices of the type (p + 
p\ )/2, generally leads to an improvement of the value 
of E (p, L) but worsens the convergence of the first ex- 
cited state Ei(p,L) during the FSM iterations. Still the 
behaviour of numerical convergence is analogous to that 
found for targeting with a single state only. Table || shows 
the DMRG calculation, first the ISM starting from L = 
12 and then the FSM iterations for the branching- fusing 



6 Finite-size scaling of the relaxation time 

We want to test the accuracy of the DMRG method ap- 
plied to the critical region of reaction-diffusion processes. 
The main interest of this study is not so much to obtain 
an extremely precise value for some exponent through a 
huge computational effort, but rather to get some insight 
into the generic behaviour of the dmrg in this kind of 
application. 



6.1 Diffusion-annihilation 

We analyse the diffusion-annihilation model defined in 
(Oh. The ground state is twofold degenerate Eq(p,L) — 
E~i{p, L) = 0, even for L finite, since the reaction AA — > 00 
reduces the number of particles and thus there are two 
stationary states: the empty lattice 1 000 ... 0) and a state 
obtained from the following combination of one-particle 
states: \A%% ...) + \%A% ...) + .. . + 1000 ...A). 

The gap is known exactly (see |ll[ and references therein) 



r(p,L)=E 2 (p,L) = 2p 1-cos 



(33) 



In the thermodynamic limit L — > oo one has r(p, L) 
1/L 2 , i.e. a dynamical exponent 8 = 2. 

Figure [l] shows a plot of the relative error: 



ML) ■■= 



E™* G (p,L)-E 2 (p,L) 



E 2 (P,L) 



(34) 



The calculation is done following the method described 
at the end of section y: The finite system method sweeps 
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Fig. 1. Relative error on the gap r(p,L) of diffusion- 
annihilation with p = 0.5 as a function of the system length 
L. For L = 16, 24, 32 and 40 (marked by the dotted vertical 
lines) we report the two values of the errors before and after 
the application of the finite system method (see text). 



Table 5. Finite-size estimates of the exponent 8 for the model 
(ll2j) as a function of the system size L. 



L 


9(L) 


10 


1.8219 


18 


1.8959 


26 


1.9265 


34 


1.9432 


42 


1.9537 


50 


1.9609 


58 


1.9662 



are applied for chains of lengths L — 16, 24, 32 and 40 
during the same run. For these sizes, in the figure two 
values of A%(L) are given, which are obtained at the be- 
ginning and at the end of the application of the FSM. We 
always observed a strong decrease of the error during the 
FSM iterations. For example, for L — 24 the point A in 
figure corresponds to a relative error equal to 10 -2 and 
was reached before the application of the FSM , while after 
the FSM iterations, we are at B, where the error drops to 
10~ 9 . Notice that a single DMRG step with the ISM from 
L = 24 (B) to L = 26 (C) corresponds to a rather big 
increase of the error from 10~ 9 to 10" 4 . 

We conclude that for critical systems, as in the case at 
hand, the use of the FSM leads to a remarkable increase 
in the accuracy, even for chains of small lengths. That is 
in general not the case for DMRG applied to symmetric 
problems. 

Furthermore, the additional numerical precision gained 
from the FSM with respect to the ISM used alone is needed 
if precise information on the critical parameters is desired. 
To illustrate this, consider the numerical calculation of the 
dynamical exponent 9. Finite-size estimates can be ob- 
tained from 9(L) = - ln(P(L + 2)/T(L))/ln((Z + 2)/L) 
and are listed in table H. While the sequence clearly con- 



verges toward the exact value 0(oo) = 2, it is also apparent 
that the finite-size data themselves are, even for the rela- 
tively large sizes used here, still quite far from the L — » 00 
limit. In conclusion, it is not possible to simply take some 
value 9(L) calculated on a large lattice to be a reliable 
estimate for the exponent 9. Rather, a careful extrapo- 
lation of the data towards their L — > 00 limit must be 
performed, as will be further discussed below. However, 
finite-lattice extrapolations are only possible if the data 
are accurate to at least 6 or 7 digits. Since even in the 
simple model (n2h finite-size effects are considerable, this 
should be expected to hold to an even larger extent in 
more complicated systems. 

The data in figure [j] were obtained from the density 

matrix p 2 , defined in (plf), with m = 32 states kept. At 
least for the examples studied in this paper, we found 
that the numerical accuracy cannot be further improved 
by increasing m. This is so since the eigenvalues u)i of the 
density matrix vanish so rapidly with increasing i that the 
truncation error is of the same order of magnitude as the 
numerical errors from the other parts of the calculation. 

For systems of sizes larger than L w 50, our results 
become numerically unstable. This might be due to the 
fact that the gap vanishes rather rapidly in this example 
(9 — 2). On the other hand, it is well-known|19| that even 
the symmetric DMRG cannot describe critical systems in 
the L — ► 00 limit, where the gap becomes too small. 



6.2 Branching-fusing 

The branching-fusing model defined by the rates (O) is 
critical at p = p c , the value of which is not known exactly. 
This critical point can be extracted from finite-size data 
of the gap, using the following scaling form, valid in the 
vicinity of p c 



r(p,L)^L- G(\p- Pc \L 1 ^ 



(35) 



where G(x) is a scaling function. Form a comparison of the 
gaps of three consecutive sizes, say L — 2,L and L + 2 one 
identifies p c (L) as the value of p for which the equation 

log [r(p, l + 2)/r(p, l)} _ log [r( P , L)/r( P , l - 2)} 



log[L/(L + 2)} 



log[(L 



2)/L] 
= 9(L) (36) 



is satisfied. In addition, an estimate 9(L) for the exponent 
9 is obtained [f23| . 

The gaps r(p, L) have been calculated for chains of 
lengths up to L rs 60, usually with to = 32 states kept. 
As before, we find that a larger value of to does not im- 
prove the results, since the density matrix eigenvalues u>i 
fall off rapidly with i and typical truncation errors for 
to = 32 are of the order e w 10~ 15 . On the other hand, 
the use of more than 64-bit arithmetic does improve preci- 
sion, which means that arithmetic precision is the limiting 
factor here. Since for an accurate extrapolation, we need 
precise finite-lattice data we merely considered chains of 
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lengths L < 30. These data are already sufficient to our 
goal to investigate the generic behaviour of the dmrg ap- 
plied to the calculation of critical parameters. 

In table pi we display the finite-size data for the critical 
point p c and several critical exponents which we are going 
to analyse. 

We begin by estimating p c . By looking at the data, we 
see that the differences p c (L + 2) — p c {L) decrease only 
slowly with increasing L. This already implies that the 
limit p c = p c (oo) must be quite far away from the finite- 
lattice data. In fact, this situation had already been en- 
countered before for periodic boundary conditions, when 
studying the finite-size scaling of Reggeon field theory, 
which is also in the universality class of directed percola- 
tion. There, it was shown [M that the L — > oo extrapo- 
lation may be reliably carried out with the BST algorithm 
pq,E6[. This algorithm transforms a 'logarithmically con- 
verging' (e.g. [J10[) sequence into another one which is ex- 
pected to converge faster. It involves a free parameter to 
which roughly measures the effective leading correction 
exponent and is chosen to achieve optimum convergence. 
Comparison of the behaviour of the BST algorithm and 
several alternative extrapolation schemes applied to finite- 
size data [ p6[p7| has shown that in the generic case when 
the correction exponent has a non-integer value, the BST 
algorithm is the most reliable of the schemes presently 
available. 

The extrapolation procedure is illustrated in table [fl. 
The first column gives the values of p c {L) which solve 
eq. (pq), taken from table @. The subsequent columns 
present the convergence-accelerated sequences as found 
from the BST transformation. For the chosen value of u>, 
stability is found and a conservative estimate of the loca- 
tion of the critical point is p c = 0.84036(1). Notice that 
the final value of p c is indeed quite far from the finite- 
lattice data, however, the fact that the given sequence 
can be made to converge well indicates that the original 
set of dmrg data is fairly accurate. The possibility of es- 
timating p c precisely mainly depends on the length of the 
sequence available for extrapolation and not so much on 
the distance of the finite-size data from their L — ► oo limit. 



In the same way, the data obtained from (|36|) for the 
dynamical exponent 9 can be analysed. Its determina- 
tion is independent of the final estimate of p c . We do 
not present the details, but simply quote our result 9 = 
1.580(1). First of all, it is satisfactory to see the good 
agreement with the value quoted in table |l|. Second, we 
point out that even for L = 28, the finite-lattice estimate 
is still some 20% away from its limit value. For compar- 
ison, we recall that for periodic boundary conditions (in 
the Reggeon field theory and with an accuracy for 9(oo) 
comparable to the situation at hand) for L = 14, the cor- 
responding difference is of the order of 3% |24|] . This is a 
consequence of the free boundary conditions usually em- 
ployed with the dmrg. It is remarkable that in spite of 
the extra difficulty presented by the free boundary condi- 
tions, one is still capable to determine 9 so precisely. In 
order to improve the precision, one would have to perform 



the DMRG with enhanced numerical accuracy in order to 
generate longer sequences. 



Differentiating numerically cq. (35), one finds 



dr(p,L) 
dp 



= L- e+1 '^ G' ( 



-P^L 1 '^ 



P~p. 



(37) 



which allows to estimate the exponent (, — 9— l/Vj.- Here, 
the numerical derivatives were calculated at the values of 
p c {L) given by the solutions of eq. (|36[). From numerical 
extrapolation we find Q = 0.66(2), and with the value of 
9 found above we find for the (spatial) correlation length 
exponent i/± = 1.08(2). This is also in good agreement 
with the value of table [l]. It is somewhat less accurate, 
however, than the estimated value for 9, since its determi- 
nation involves the calculation of numerical derivatives. 
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Fig. 2. Scaling plot of the particle density L*' v± n(l) as a 
function of the scaled distance l/L at the estimated critical 
point p c — 0.84036 and for an injection rate of p' = 0.3. 



7 Density profiles in the branching - fusing 
model 

7.1 Finite-size scaling of the density profiles 

Besides the calculation of gaps, the dmrg allows investi- 
gating density profiles in the steady state. If h(l) is the 
density operator at position I, the density profile is nat- 
urally calculated from the ground states \s) = \ipQ ) an d 



M 



q |, normalised as (s\s) 



1, as 



n(!) = (1>iPW)\4 r) ) 



(38) 



However, for the models (h2|Jl3|), this will be non- vanishing 
only if particles are injected at the boundaries. Therefore, 
we add the boundary reaction (|14J) to the branching-fusing 
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model and analyse the behaviour of the system as a func- 
tion of the injection rate p' . 

Adapting the finite-size scaling arguments formulated 
originally by Fisher and de Gennes pa ] for equilibrium 
profiles, the density of particles at the critical point p c 
should follow a scaling form 



i(l) = rM v ±F(l/L) 



(39) 



where < / < L denotes the position along the system, 
F(x) is a scaling function and [3 the order parameter crit- 
ical exponent. According to (|39] ) the quantity L /3 ^ u± n(l) 
depends on I only through the scaling variable l/L. 

Fig. g shows a scaling plot of L^^ u± n(l) at p = p c 
as determined above and for particle injections rate p' = 
0.3, using the value (3/v±_ = 0.25208, as quoted in J17J. 
We see that the densities for systems of lengths L — 
20, 24, 28 ... 56 tend to collapse nicely onto a single curve. 
Nevertheless, notice that finite-size corrections are quite 
large, in particular, they are larger than encountered for 
example in Ising model calculations. 

We now try to find p c and (3/v± from our finite-lattice 
data. Consider the point I = L/2 in the middle of the 
chain. The central density n(L/2;p) should obey in the 
vicinity of the critical point the scaling law 



i(L/2;p) = L-V^H (\p - p c \L 1/v± 



(40) 



where H (x) is a scaling function. In analogy to the analysis 
of the gap, we can compute finite-size estimates for the 
critical point p c and the exponent f3/v± by replacing in 
eq. (B6) r by n and 9 by /3/v±, with the results listed in 
table^, for p' — 0.3. 

Table H shows the BST extrapolation table for the val- 
ues of p c (L), for which we used strips of widths up to L = 
40. The extrapolation yields p c = 0.8406(3). While this 
agrees quite well with the earlier estimate p c = 0.84036 
found from the gap r, the apparent difference between 
the two results gives an a posteriori assessment of the re- 
liability of the extrapolation procedure. 

In addition, we observe that the sequence of values 
for p c obtained from r increases with increasing L, while 
the sequence obtained from n decreases. Although the raw 
data in table M are much farther away from the L — > oo 
limit of p c than those in table @, the extrapolation in the 
former case can be carried out to a higher degree of preci- 
sion than in the latter case. This illustrates once more the 
importance of a careful finite-lattice extrapolation proce- 
dure when trying to extract precise parameter values from 
lattice calculations. 

In a similar fashion as done for the exponent 8, we have 
estimated the ratio /3/v± = 0.24(1). These results, both 
for p c and for the exponents, are less accurate than those 
obtained from the finite-size scaling analysis of the gap, in- 
dicating that the numerical accuracy of eigenvectors of the 
non-symmetric Hamiltonian H for the branching-fusing 
model is inferior to that of its eigenvalues. 



7.2 The limit p ' ->• 



We now discuss the consequences of varying the boundary 
injection rate p' , in particular the limit p' — > 0. 

Fig. shows the particle density n(l) as function of the 
scaled variable l/L for an injection rate p' = 0.03, L = 16, 
24, 32, 40 and 48 at the estimated critical point p = p c . 
The behaviour is now completely different with respect to 
the profiles with p' = 0.3. Apart from the two boundary 
sites (corresponding to I = and I = L), where the den- 
sity of particles is weakly dependent on the system size L, 
we see that starting from the edges the density increases 
towards the bulk. For small systems the maximum of the 
density is found in the middle of the chain, while for sys- 
tems large enough two maxima are located at a distance 
(max w 7 lattice spacings from the edges (we find that l max 
is independent of L for sufficiently long chains). 

This phenomenon is rather counterintuitive since one 
would expect that the density of particles as a function of 
I should decay monotonically starting from the edges and 
moving towards the center, as seen in figure @ for p' = 0.3. 
The effect observed here actually has a counterpart in the 
magnetization profiles of equilibrium spin systems in the 
presence of a weak surface magnetic field [g9| . This effect 
in turn was found and explained by appealing to the uni- 
versal short-time critical dynamics involving the so-called 
slip exponent |3Q]. Non-monotonous profiles of this kind 
are a true fluctuation effect and cannot be explained on 
the level of a mean-field approximation, see ps| , p0| , pl| , p2| 
and appendix |C| It is known, for instance for the Ising 
model at its critical point, that a weak surface field hi 
induces a macroscopic scale length l\ such that up to a 
distance I w l\ from the surface, the magnetization is in- 
creasing with I. It starts to decrease again as I > 1%, with 
the asymptotic behaviour ~ l~^l v for large distances (in 
the Ising model, z/j_ = v\\ = v) fl33"[. 
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Fig. 3. Particle density n(l) at the estimated critical point 
Pc — 0.84036 for an injection rate p — 0.03 plotted as function 
of l/L. 
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To describe non-monotonic densities of particles, it is 
necessary to modify the scaling form ( |39J ) to include a term 
depending on the boundary reaction rate p' , following the 
standard theory of boundary critical phenomena at the 
ordinary transition, see [ p9|3l| . The rate p' will enter into 
the particle density in the form of a scaling variable A = 
p'l Xl , where x\ is some scaling dimension. One expects the 
scaling form 



n^p^^r^^fil/L^p'l* 1 ) 



(41) 



Following the same analysis as presented in J2!||51|, one 
finds that x\ = j3i/v_\_, with (3\ the (ordinary) order pa- 
rameter surface exponent. This exponent has been calcu- 
lated for directed percolation, see table |l|, yielding x\ ~ 
0.669. 

For large A, one should recover from (^TJ) the scaling 
form given in eq. ( p9| ) , thus 



lim T{l/L,\) 



F(l/L) 



(42) 



More interesting, for the present case, is the other limit 
where A — > 0. In absence of particle injection (p' = 0), the 
particle density n vanishes, which implies T[x, A = 0) = 0. 
At small p', it is natural to assume that the density varies 
linearly with p', leading to the condition 



F{l/L,\)~\ 



(43) 



valid in the range, say, < A < Ao- Notice that Ao defines 
a length scale l\ = (Xo/p') 1 ^ 1 ■ Therefore, in the L <C ^i 
limit, one finds 



n{l,L) = l p T(l/L) 



(44) 



obtained by inserting the limiting value (|43|) in eq. ( [4l|) 
and where p' = x\ — [3/v± = (Pi — /3)/V_i_, in complete 
analogy with |2|] . 

From the numerical values quoted in table |lj, we find 
p' ~ 0.416, the positivity of which explains the increase 
of the profiles with I close to the boundary. Actually, the 
particle density with p' = 0.03 shown in Fig , ^ corresponds 
neither to the scaling regime ( p9[ ) nor to (f44|), but to an 
intermediate situation. 

In figure ^, we display the critical particle density for 
L = 16, 20, ... 40 but with a decreased injection rate 
p' = 0.002. Profiles are monotonously increasing from the 
edges and there is no sign of an inflection of the curves 
as seen in Fig. ||. The inflection point should be found at 
a distance from the edges l max (p' — 0.002) given by the 
relation: 



ima*(p' = 0.002) = l max (p' = 0.03) 



0.03 
0.002 



l/zi 



400 



in units of the lattice constant, and where we have used 
lmax(p' = 0.03) ~ 7 as estimated from the profiles of Fig. 
H. Thus, we expect that for p' = 0.002, at about 400 lattice 
spacings from the surface the density profiles should show 
a maximum and then start decreasing again. This effect 



should be noticeable only for systems of sizes L > 2/ max ps 
800, well beyond the range of sizes shown in figure 0. 

The inset of figure Q shows the scaled profiles L~ p n(l) 
plotted as a function of the scaled distance l/L (we have 
used the expected value of p' = 0.416 for directed percola- 
tion and removed the edge sites I = and I = L from the 
plot). If the scaling assumptions leading to eq. ( |44| ) were 
correct, all the profiles for different values of L should col- 
lapse into a single curve. Indeed, we find a trend towards 
a data collapse, although finite-size corrections are rather 
strong, especially in the middle of the chain, while they 
are weaker at the edges. It is remarkable, that the scal- 
ing law (|44|) , whose derivation involves the scaling close to 
the boundary, should be valid for the entire finite system, 
at least for the lattice sizes under consideration here. On 
the basis of these observations, one may conclude that the 
DMRG data for the profiles agree with the results of the 
scaling theory also in the weak injection rate regime. 



7.3 Profiles from matrix elements 

So far, the calculations of the density profiles used the 
straightforward form ( p8[ ) and then tried to perform the 
limit p' — > numerically, which in principle should be 
taken only after the L — > oo limit has been carried out. 
Alternatively, we may rely on an analogy with the calcu- 
lation of the order parameter of equilibrium spin systems 
which avoids this cumbersome double limit ]34| , see also 

P0[ . If \ipi ) and (ip\ '\ are the first excited eigenstates of 
H, consider fB5| 



Mi 



N(l) := (^ln(Ol^r) 



«.W\ 



(45) 



where n(l) is the density operator at position I and where 
we have simply set p' — 0. Although N(l) is not directly 



0.08 




0.06 



0.04 



Fig. 4. Density of particles at p c and for p' = 0.002. The inset 
shows the collapse of the scaled densities onto a single curve 
as expected from eq. (M). 
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related to the more "physical" density n(l), it offers the 
distinct advantage that it is non- vanishing even for p' = 
and furthermore displays the same scaling behaviour as 
expected for n(l). Since in the percolating phase, the first 
gap r is exponentially small for L large, one might expect 
to find for N(l) the same scaling behaviour as for n(l). 
Indeed, as we shall see, the following finite-size scaling 
behaviour holds for N(l) at p = p c : N(l) ~ L~ /3 /" ± for 
sites deep in the bulk, i.e. l/L ~ 1/2 and N{1) ~ L~^^ v± 
for boundary sites, i.e. I = or I = L; with the same 
exponent values as found before. 

Figure g] shows a scaling plot of L@/ U± N{1) as a func- 
tion of l/L. where we used the value of (3/v± = 0.25208 
from table [y. The data collapse in a satisfactory way in the 
bulk, but rather poorly at the surface. This is consistent 
with our expectation of two different scaling regimes. It is 
also quite remarkable that the numerical data for I = L/4 
and I = 3L/4 show an almost perfect collapse. 

We analyzed the numerical data at p — p c forming 
the finite-size estimates from two successive system sizes 
as done in eq. (p6|), using the values of N(l) at / = L/2 
and at the surface 1 = 1. The data are listed in table |[ 
Extrapolation for L — ► oo with the bst method yields 
P/v x = 0.249(3) and @i/v± = 0.667(2), both in good 
agreement with the results quoted in table GL Apparently 
the exponents extrapolated from N(l) are more accurate 
than those obtained from the density profiles and allow a 
direct determination of surface exponents, allowing a com- 
parison with field-theory methods |3§] . Again, the data for 
N(l) might also be used to find yet another estimate of 
p c , but we have not carried out this calculation. 



8 Conclusions 

In this paper we have investigated the properties of two 
reaction- diffusion models at their critical point, by means 



1.2 




0.4 



0.2 



L = 36 



0.4 0.6 

l/L 



Fig. 5. Plot of the scaled quantities L"' v± N(l) as function 
of the scaled distance l/L for L — 16, ... at the critical point 



P=Pc 



of DMRG techniques. Our study demonstrates that the 
dmrg is capable to treat rather well out-of-equilibrium 
systems and provides a method for calculations of critical 
exponents, alternative to other general methods such as 
Monte Carlo simulations or series expansion techniques. 
The dmrg offers the advantages that (i) there is no criti- 
cal slowing-down, (ii) the method does not require random 
numbers and (iii) there is no need to make any assump- 
tions about the basic state around which one may expand. 
Our findings may be summarised as follows. 

(1) In the examples considered, the symmetric density ma- 
trix pM produced the most accurate results. It is essential 
to have a reliable method to diagonalize a sparse non- 
symmetric matrix. The non-symmetric Lanczos and the 
Arnoldi algorithms used here were found to produce re- 
sults of comparable accuracy. 

(2) It was enough to keep to — 32 states. Improvements 
in numerical accuracy of the dmrg cannot be achieved 
by increasing to, but rather by enhancing the number of 
digits kept beyond the conventional 64-bit arithmetics. 

(3) In order to analyse the critical region, it is essential 
to use the finite system method of the dmrg. The infinite 
system method alone is not sufficient. 

(4) Finite-size data for critical non-equilibrium systems 
are strongly affected by finite-size corrections. Because 
usually the dmrg works best for free boundary conditions, 
finite-size corrections should be expected to be substan- 
tially larger than in calculations done with standard diag- 
onalization techniques, which usually work with periodic 
boundary conditions. However, the trade-off is that from 
a single DMRG calculation one may get both bulk and sur- 
face exponents. 

(5) In studies of the scaling of profiles, the consideration of 
matrix elements of the type of ( |45| ) may provide a better 
computational tool than the physically more straightfor- 
ward profiles (|38|). 

(6) Precise estimates of the location of the critical point 
and of critical parameters can only be obtained through 
a finite-sequence extrapolation technique. That requires a 
long sequence of precise lattice data for several distinct 
sizes, rather than data from a single large lattice. 

(7) The extrapolation can be carried out reliably by the 
bst algorithm. The results for the critical exponents are in 
agreement with the results for series expansion and Monte 
Carlo simulation, see table [j]. At present, the DMRG do not 
yet achieve the same precision as these. The precision of 
the method can be improved by generating longer chains. 
To overcome the numerical instabilities encountered, this 
will require to go beyond the usual 64-bit arithmetics. 

All in all, we conclude that the DMRG has turned out 
to be a reliable general-purpose method, also for study- 
ing non-equilibrium critical phenomena, although it works 
best out of criticality. 
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A Relation with site and bond percolation 

We recall the relationship of the reaction-diffusion pro- 
cesses defined in w\-iw with directed percolation. The lat- 
ter is usually defined in terms of local probabilities giving 
the state of a site at time t + 1 in terms of its two parent 
sites at time t 







PsPb 



PsPb (2- Pb ) (46) 



where p s and p b describe site and bond percolation, re- 
spectively. Here • marks an occupied and o an empty 

site. Furthermore, = 1 — and so on. From these 

o • 

rules, a time evolution operator S, acting according to 
\P(t + 1)) = S\P(t)), is constructed, which is related to 
the quantum Hamiltonian H = 1 — S. 

In order to link the percolation parameters p s , p b with 
the reaction and diffusion rates a, (3, 7, 8 and D, we ob- 
serve that the evolution of a site is according to ( fig ) only 
dependent on its parentsand independent of the state of 
its neighbours at the same time t. Thus 



• o 



• o 

• o 



• o 



• o 
o • 



• o 



(47) 



by summing over the possible states of the right or left 
nearest neighbours. Similar relations hold for the other 
rates. From this, the following conditions for a mapping 
of the process (|l|-||) onto directed percolation hold 

P = l-8=p s pb ; D = ; 2a + 7 = 1 - p s p b (2 - p b ) 

_ (48) 

It is easy to check that the model (|l3|) cannot be exactly 
mapped onto directed percolation this way, although it 
still is in the same universality class. For example, a map- 
ping onto site percolation (p b = 1) requires that 5 — 
2a + 7, while bond percolation (j> s = 1) is achieved for 
5 2 = 2a + 7. 



B On the density matrix 

To justify the use of the density matrix tr |"0o) {^oIj onc 

writes the targeted state in a product state of states | ai , ii) = 
\i) of the "system" one wants to describe and \j r , (3 r ) = \j) 
of the "environment" . From the construction of the dmrg, 
each of the bases has mn elements, where n is the number 
of states per site. Then 



where the sum over a runs over m orthogonal states in 
the system basis, \a) = ^. u a i\i). One now demands that 
the approximation is optimal by minimising 



liM = X>(»,i)|»)|j) 

One now approximates \ipo) by \ipo), with 
IV>o) = '^''Poioi, j)\a)\j) 



(49) 



(50) 



||Vo)-|^o)|| 2 = l-2^^ 2 (a,j). 



(51) 



Making this stationary with respect to il>a(a,j) gives 

Si ^o(*7 j)u a i = ipo(<x,j) and thus the global minimum 
of 

II |Vto) " IV'o) f = 1 - Y, U °"P^ *')««' (52) 

Oiii' 

has to be found. Here one introduces the reduced density 
matrix 

P(M') = ^2ipo{i,j)ipo{i',j)- (53) 

i 

The stationary points of (p2) are obviously given by set- 
ting \a) equal to the eigenvectors of the reduced density 
matrix (Ritz condition). The stationary values of ( |52| ) are 
then given by 1 — ^Z i= x \i, where Ai are the eigenvalues 
of the density matrix (0 < Aj < 1, J^ Aj = 1). The global 
minimum is obtained by choosing the m eigenvectors with 
the largest associated eigenvalues. 

In the non-symmetric case, repeating the same argu- 
ment, one finds that the minimization of: 



-,(0\ 



',(' 



\m')-m'}\r + \M')-\t 



M\ 



„7.( r )\l|2 



(54) 



yields as optimal basis set the eigenvectors of the density 
matrix (BQ). 



C Density profiles from kinetic equations 

We show that for branching- fusing (131) , the mean field 
steady-state density profile for the semi-infinite system 
x > with a prescribed boundary density is a monotonous 
function of the distance x from the boundary. If a (a;, t) is 
the mean particle density, the kinetic equation is 



Da" + 4(p - 1/2) a - 2a 2 



(55) 



In the steady state, a = 0. We write a(x) = aoo^(^/^j_), 



where a r 



2p — 1 ~ £ , is the bulk density and £± 



yjD/(p — 1/2) the spatial correlation length. The profile 
is 



3 / V(2y + l)/3 + tanh(y) \ _1 
nV) 2 \l + ^/(2 i p + l)/3 tanh(y)i 2 



(56) 



where <po = a(0)/a oo is related to the boundary den- 
sity. Evidently, tp(y) — > 1 as y — > 00 monotonously. The 
above result was derived for a semi-infinite system with 
p > PcMF — 1/2. For a finite system at p = 1/2, the 
p-dependence of the correlation length £± is traded for a 
size-dependence. 
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Table 3. Gaps r = E\(p,L) for the branching-fusing model eq. (|13|) for p = 0.8403578 and with m = 32 states kept as found 
from the iterations of the fsm using the density matrices p' 3 ' and p' 1 '. Li and L r indicate the lengths of the left and right 
parts, which vary during the application of the finite system method. The last column was calculated using 30 digit precision 
arithmetic, yielding a result free of diagonalization inaccuracies. 



U 


l^r 


r from p l3J 


r from p |1J 


r from p' 1 ', high prec. 


6 


6 


0.0211852795111 


0.0211852795111 


0.0211852795111 


7 


7 


0.0173940529006 


0.0173940538620 


0.0173940538302 


8 


8 


0.0146003381516 


0.0146003960454 


0.0146003961355 


9 


7 


0.0146004093819 


0.0146003960889 


0.0146003961161 


10 


6 


0.0146006193716 


0.0146003961644 


0.0146003961865 


9 


7 


0.0146004129822 


0.0146003961642 


0.0146003961866 


8 


8 


0.0146004043568 


0.0146003962502 


0.0146003962379 


7 


9 


0.0146004373237 


0.0146003963897 


0.0146003962378 


6 


10 


0.0145993594105 


0.0146003962420 


0.0146003962377 


7 


9 


0.0146004306753 


0.0146003961645 


0.0146003962378 


8 


8 


0.0146004003523 


0.0146003961420 


0.0146003962379 


9 


9 


0.0124729781067 


0.0124729862108 


0.0124729862559 


10 


8 


0.0124730375718 


0.0124729861363 


0.0124729862333 


11 


7 


0.0124731467646 


0.0124729861847 


0.0124729862262 


12 


6 


0.0124732937813 


0.0124729861234 


0.0124729862244 


11 


7 


0.0124731477916 


0.0124729861816 


0.0124729862251 


10 


8 


0.0124730400858 


0.0124729861765 


0.0124729862258 


9 


9 


0.0124730036975 


0.0124729869859 


0.0124729861961 



Table 4. Gaps _T = E2(p, L) for the diffusion- annihilation model eq. (fl2|) for p — 1 and with m = 32 states kept as found from 
the iterations of the fsm using the density matrices p' 3 ' and p"'. Li and L r indicate the lengths of the left and right parts, 
which vary during the application of the finite system method. The last column gives the exact result. For total length 18, only 
the three results for equal lengths are given. 



L, 


L r 


r from p l3J 


r from p |1J 


exact 


6 


6 


0.0581163650912 


0.0581163651488 


0.0581163651479 


7 


7 


0.0437254502992 


0.0437047728086 


0.0437047985324 


8 


8 


0.0340666763278 


0.0340537551881 


0.0340538006322 


9 


7 


0.0340537973881 


0.0340537521644 




10 


6 


0.0340538348398 


0.0340537779114 




9 


7 


0.0340437188607 


0.0340537779110 




8 


8 


0.0340537940475 


0.0340538006303 


0.0340538006322 


7 


9 


0.0340537902152 


0.0340538006306 




6 


10 


0.0340537976329 


0.0340538006308 




7 


9 


0.0340537982237 


0.0340538006312 




8 


8 


0.0340537985217 


0.0340538006323 


0.0340538006322 


9 


9 


0.0272774311113 


0.0272774426890 


0.0272773931946 


9 


9 


0.0272773902380 


0.0272773931864 


0.0272773931946 


9 


9 


0.0272773945078 


0.0272773931999 


0.0272773931946 
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Table 6. Finite-size estimates of critical point y r and of various exponents for the branching-fusing model, obtained from the 
gap r, from the density profile n(l) defined in (|38J ) and from the scaUng of N(l) as defined in (|45]). The last row shows the 
L — » oo limit obtained from bst extrapolation. The numbers in brackets give the estimated uncertainties in the last digit. 





from r 


from 


n(l) 


from N(l) 


L 


Pc(L) 


6{L) 


cw 


Pc{L) 


P/v±(L) 


0/v±{L) 


Pi/v±{L) 


10 


0.815486295 


0.830071389 


0.177917024 










12 


0.822241704 


0.923515450 


0.248868030 






0.211498060 


0.524156106 


14 


0.826556808 


0.996672190 


0.303005003 






0.214641534 


0.540022433 


16 


0.829477408 


1.055258740 


0.345372694 


0.844595690 


0.174469664 


0.217449273 


0.552833181 


18 


0.831547147 


1.103159519 


0.379339662 






0.219928645 


0.563377349 


20 


0.833068754 


1.143030157 


0.406998719 


0.843578941 


0.183190533 


0.222113618 


0.572197559 


22 


0.834221223 


1.176727177 


0.430122954 






0.224042959 


0.579677788 


24 


0.835115836 


1.205580740 


0.449620922 


0.842813911 


0.191175959 


0.225753163 


0.586096418 


26 


0.835824726 


1.230565614 


0.466327849 






0.227274607 


0.591673398 


28 


0.836396350 


1.252411806 


0.480376056 


0.842276687 


0.197772094 


0.228648854 


0.596559909 


30 












0.229893751 


0.600879333 


32 








0.841894149 


0.203168874 


0.230988495 


0.604833883 


34 












0.231983300 


0.608060181 


36 








0.841617080 


0.207584161 


0.232890344 


0.611490572 


40 








0.841415905 


0.211176020 






oo 


0.84036(1) 


1.580(1) 


0.66(2) 


0.8406(3) 


0.24(1) 


0.249(3) 


0.667(2) 



Table 7. BST extrapolation table for the critical point p c . The first column are the "raw" data obtained from r (we used L = 16, 
18, 20, . . . 28) and the following columns are obtained by repeated application of the BST transformation, with uo = 1.8475. 



0.829477408 


0.840171029 


0.840328212 


0.840352241 


0.840366202 


0.840357906 


0.840357626 


0.831547147 


0.840223500 


0.840337691 


0.840357304 


0.840358276 


0.840357764 




0.833068754 


0.840258488 


0.840344750 


0.840357765 


0.840357763 






0.834221223 


0.840282908 


0.840349176 


0.840357764 








0.835115837 


0.840300341 


0.840351951 










0.835824726 


0.840313022 












0.836396350 















Table 8. BST extrapolation table for the critical p c from the scaling of the density of particles n(L/2) in the middle of the 
system, where sizes from L = 16 up to L — 40 were used, with u> — 2.316. 



0.844595690 


0.842080802 


0.840373866 


0.840653412 


0.840568730 


0.840563673 


0.840569610 


0.843578941 


0.841361648 


0.840593970 


0.840558615 


0.840564091 


0.840556196 




0.842813911 


0.841027230 


0.840571147 


0.840563201 


0.840572164 






0.842276687 


0.840840423 


0.840566363 


0.840610429 








0.841894149 


0.840734840 


0.840594342 










0.841617080 


0.840683057 












0.841414666 















